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We give a pedagogical review of the application of field theoretic and path integral 
methods to obtaining perturbative solutions of stochastic differential equations. 



There are many applications of stochastic differential equations (SDE) for mathematical 
modeling. In the realm of neuroscience, SDEs are utilized to model stochastic phenomena 
that range in scale from molecular transport in neurons, to neuronal firing, to networks of 
coupled neurons, to even cognitive phenomena such as decision problems pQ. In many appli- 
cations, what is often desired is the ability to obtain closed form solutions or approximations 
of quantities such as the moments of stochastic processes. However, generally these SDEs 
are nonlinear and difficult to solve. For example one often encounters equations of the form 



where r](t) represents some noise process, in the simplest case a white noise process where 
(i](t)) = and (rj(t)r](t')) = 5(t — t'). Often of interest are the moments of x(t) or the 
probability density function p(x,t). Traditional methods use Langevin or Fokker-Planck 
approaches to compute these quantities, which can still be difficult and unwieldy to apply 
perturbation theory [2H3] . Here, we will show how methods developed in nonequilibrium 
statistical mechanics using path or functional integrals can be applied to solve SDEs. 

While these methods have been recently applied at the level of networks, the methods are 
applicable to more general stochastic processes [T3HT8] . Path integral methods provide 
a convenient tool to compute quantities such as moments and transition probabilities per- 
turbatively. They also make renormalization group methods available when perturbation 
theory breaks down. 

Although Wiener introduced path integrals to study stochastic processes, these methods 
are not commonly used nor familiar to much of the neuroscience or applied mathematics 
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community. There are many textbooks on path integrals but most are geared towards 
quantum field theory or statistical mechanics [T3H2T] . Here we give a pedagogical review of 
these methods specifically applied to SDEs. In particular, we review the response function 
method j22j E3], which is particularly convenient to compute desired quantities such as 
moments. 

The goal of this review is to present methods to compute actual quantities. Thus, mathe- 
matical rigor will be dispensed for convenience. This review will be elementary. In Section [IT], 
we cover moment generating functionals, which expand the definition of generating func- 
tions to cover distributions of functions, such as the trajectory x(t) of a stochastic process. 



We continue in Section III by constructing functional integrals appropriate for the study 



of SDEs, using the Ornstein-Uhlenbeck process as an example. Section IV introduces the 
concept of Feynman diagrams as a tool for carrying our perturbative expansions and in- 
troduces the "loop expansion", a tool for constructing semiclassical approximations. The 
following section [V] provides the connection between SDEs and equations for the density 
p(x, t) such as Fokker-Planck equations. Finally, we end the paper by pointing the reader 
towards important entries to the literature. 



II. MOMENT GENERATING FUNCTIONALS 



The strategy of path integral methods is to derive a generating function or functional 
for the moments and response functions for SDEs. The generating functional will be an 
infinite dimensional generalization for the familiar generating function for a single random 
variable. In this section we review moment generating functions and show how they can be 
generalized to functional distributions. 

Consider a probability density function (PDF) P(x) for a single real variable x. The 
moments of the PDF are given by 

<*"> = /*w 

and can be obtained directly by taking derivatives of the generating function 

Z(X) = (e Xx ) = [ e Xx P(x)dx 
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with 



(x») 
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Note that in explicitly including Z[0] we are allowing for the possibility that P(x) is not 
normalized. This freedom will be convenient especially when we apply perturbation theory. 

(x-a) 2 

For example, the generating function for a Gaussian PDF, P(x) oc e ^ 2 , is 

ix — a) , * 

e--^ +Xx dx (1) 

■oo 

The integral can be computed by completing the square so that the exponent of the integrand 
can be written as a perfect square 

\2 



[x — a) 
2a 2 



+ Xx = -A(x - x c y + B 



This is equivalent to shifting x by x c , which is the critical or stationary point of the exponent 

d ( {x — a) 2 



dx 
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Xx 







yielding x c = Act 2 + a. The constants are then 
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x 2 a 2 X 2 a 2 



r> _ c 

2ct 2 2ct 2 



+ Xa 



The integral in (IT]) can then be computed to obtain 



{x-Xtr 2 -a) 2 , A 2 ct 2 , , \ , > 2 " 2 
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Z(X) = / e — +Aa+— dx = z ^ e 

J — oo 

where 

/OO 2 
■oo 

is a normalization factor. The mean of x is then given by 

(x) - " - Xa+ 



A o 



d ■ ^ 2 ~ 2 



dX 



e — ' 2 



A=0 



The cumulant generating function is defined as 



W(X) = lnZ{X) 
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so that the cumulants are 



A=0 



In the Gaussian case 

W(X) = \a+ ^AV + lnZ(0) 

yielding (x)c — (x) — a and (x 2 ) c = var(x) = (x 2 ) — (x) 2 = a 2 , and (x n )c = 0, n > 2. 

The generating function can be generalized for an n-dimensional vector x = 
{x±,X2, ■ • ■ ,x n } to become a generating functional that maps the n-dimensional vector 
A = {Ai, A 2 , . . • , A„} to a real number with the form 



/n 
n 
i=i 



where Kj£ = (K~ l )j k and we use square brackets to denote a functional. This integral can 
be solved by transforming to ortho normal coordinates, which is always possible if is 
symmetric, as it can be assumed to be. Hence, let u a and v a be the ctth eigenvalues and 
orthonormal eigenvectors of K~ l respectively, i.e. 

and 

3 

Now, expand x and A in terms of the eigenvectors with 

a 

Afc = ^ d a v% 

a 

Hence 

Y X i K Jk X * = Y Y CaUpCfiVjVj = Y ^^WW = Y UaC * 
j,k j a, f3 a, a 
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Since, the Jacobian is 1 for an orthonormal transformation the generating functional is 

Z[X\ = J Y[dc a e^-^ +daCa) 



a 



= | J / dc a e-^ aci+daCa 

a 

= Z[0}e^ k ^ XjKjkXk 

where 

Z[0] = (2ndetK) n/2 
The cumulant generating functional is 

W[X] = \nZ[X] 



Moments are given by 



Xi=0 

However, since the exponent is quadratic in the components A;, only even powered moments 
are non-zero. From this we can deduce that 

\i=l I all possible pairings 

which is known as Wick's theorem. Any Gaussian moment about the mean can be obtained 
by taking the sum of all the possible ways of "contracting" two of the variables. For example 



Xd) = K ab K cd + K ad K bc + K ac K, 



bd 



In the continuum limit, a generating functional for a function x(t) on the real domain 
t G [0, T] is obtained by taking a limit of the generating functional for the vector Xj. Let 
the interval [0, T] be divided into n segments of length h so that T = nh and x(t/h) = Xi 
for t G [0, Tj. We then take the limit of n — > oo and h — > preserving T = nh. We similarly 
identify Aj — > X(t) and Kij — > K(s,t) and obtain 



Z[\] = J -pa;(t)e _ 5/ ;c ( s )^ _1 ( s '*) :E W d ^+/ A W :E W d * 
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= Z{0]e^ x( -^ K( - s ^ x ^ dsdt (2) 
where the the measure for integration 

n 

Vx(t) = lim TTrf^i 



n— too ■ 

i=0 



is over functions. Although Z[0] = lim n _^. 00 ( 27r det K) n l 2 is formally infinite, the moments 
of the distributional are well defined. The integral is called a path integral or a functional 
integral. Note that Z[X] refers to a functional that maps different "forms" of the function 
X(t) over the time domain to a real number. Defining the functional derivative to obey all 
the rules of the ordinary derivative with 

5X(s) 



SX(t) 

the moments again obey 



S(s - f) 



For example 



K(t h ,t i2 )---K(t i2s _ i: t ti J 

all possible pairings 



W'^=w\A)sm zlx]=K(h ' k) 



We can further generalize the generating functional to describe the probability distribution 
of a function ip(x) of a real vector x, instead of a single variable t with 

Z[X] = j ' Vpe~^ 'v(y)K- 1 (y,xMx)d d yd d x+f \(g)<p(x)d d x 

= z[0}e^ x{g)Ki ^ )x{£)ddyddx 

Historically, computing moments and averages of a probability density functional of a func- 
tion of more than one variable is called field theory. In general, the probability density 
functional is usually written in exponential form 



P[<P] = e 



-sir®] 
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where S[ip] is called the action and the generating functional is often written as 



where 




For example, the action given by 

S[<p] = J v{£)K- l {t,t')v{t')d d td d t' + g J ^ A {f)d d t 

is called (p A ("<p-4") theory. 

The analogy between stochastic systems and quantum theory, where path integrals are 
commonly used, is seen by transforming the time coordinates in the path integrals via 
t — 7- it (where i 2 = —1). When the field is a function of a single variable t, then this would 
be analogous to single particle quantum mechanics where the quantum amplitude can be 
expressed in terms of a path integral over a configuration variable (f)(t). When the field is 
a function of two or more variables 0(r, t), then this is analogous to quantum field theory, 
where the quantum amplitude is expressed as a path integral over the quantum field </>(f, t). 



III. APPLICATION TO SDE 

Building on the previous section, here we derive a generating functional for SDEs. Con- 
sider a Langevin equation 

^ = f(x,t)+g(x,t) V (t) 

with initial condition x(t ) = y, on the domain t G [0, T\. Equation ^ is to be interpreted 
as the Ito stochastic differential equation 

dx = f(x,t)dt + g(x,t)dB t (3) 

where dB t is a Brownian stochastic process. We will show how to generalize to other 
stochastic processes later. According to the convention for an Ito stochastic process, g(x,t) 
is non- anticipating, which means that in evaluating the integrals over time and B t , g(x, t) is 
independent of B T for r > t. The choice between Ito and Stratonovich conventions amounts 
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to a choice of the measure for the path integrals, which will be manifested in a condition on 
the linear response or "propagator" that we introduce below. 

The goal is to derive a probability density functional (PDF) and moment generating 
functional for the stochastic variable x(t). For the path integral formulation, it is more 
convenient to take x(to) = in ^ and enforce the initial condition with a source term so 
that 

dx 

— = f(x,t)+g(x,t) V (t)+y8(t-t ) (4) 

where S(-) is the point mass or Dirac delta functional. The discretized form of Q with the 
Ito interpretation for small time step h is given by 

x i+ i - Xi = fi(xi)h + gi(xi)wiVh + y5 it0 (5) 

j G {0,1,..., N}, T = Nh, 5ij is the Kronecker delta, xq = 0, and W\ is a discrete random 
variable with (wi) = and (wiWj) = S i: j. Hence, the discretized stochastic variable vector 
Xi depends on the discretized white noise process Wi and the initial condition xq. We use x 
and w without indices to denote the vectors x = (xx, . . . , x^) and w = (wq, Wi, . . . , tOzv-i)- 
Formally, the joint PDF for the vector x can be written as 

N 

P[x\w; y] = Yl 5[x i+1 - Xi - fi(xi)h - gi(xi)wiVh - y5 ifi ] 

i=0 

i.e. the probability density function is given by the point mass (Dirac delta) constrained at 
the solution of the SDE. 

Inserting the Fourier representation of the Dirac delta 

gives 

N 

Lp-i-Y^i kj(x j+1 -x j -f j (x j )h-g j (x j )w j Vh-y8j fi ) 



P[x\w;y] = ff[ 



3=0 

The PDF is now expressed in exponential form. 

For Gaussian white noise the PDF of u>j is given by 

P(v)i) = -^e- 1 ^* 
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Hence 



f n 
P[x\y) = / P[x\w;y) \\P{w j )dw j 

J ,;—n 



3=0 

= f JJ^?. e -*Ei*i(*i+i-*i-/j(*j)fc-»«j,o) f >ll ' J ( ^ H 
can be integrated by completing the square as demonstrated in the previous section to obtain 

Taking the continuum limit h — > 0, N oc such that T = iV/i gives 

P[x(t)\y,t ) = [ P5(t) e ~/[ £ W( i W~ / ( ;c W' t ^ s/<5 (*-* ))-^ 2 W s,2 ( :r W' i )] d ' 



with a newly defined complex variable iki — >■ x(t). We include the argument of x(t) as a 
reminder that this is a functional of x conditioned on two scalars y and to- The moment 
generating functional for x(t) and x(t) is then given by 



Z[J, J) = J Vx{t)Vx{t)e 



S[x,x]+f J(t)x(t)dt+f J(t)x(t)dt 



with action 

1 



S[x, x] 



x(t)(x(t) - f(x(t),t) - y5{t - t )) - ^x 2 (t)g 2 (x(t), t) 



dt (6) 



The probability density functional can be derived directly from the SDE (J3]) by considering 
the infinite dimensional Dirac delta functional and taking the path integral: 

P[x(t)\y,t Q ] = [ V V (t)5[x(t)-f(x,t)-g(x,t) V (t)-y5(t-t )}e-^ 2 (t)<it 

Vg{t)Vx{t)e~^^^~ f ^ x ^ 
j)^ e - f m(Ht)-f^,t)~yS(t-t ))+h 2 xHt)g 2 (x,t)dt 

yielding the action (|6]) [29J. Owing to the definition iki — > x(t) the integrals over x(t) are 
along the imaginary axis, which is why no explicit i appears in the action above. 

In a similar manner, we can define the path integral for more general processes than the 
Brownian motion processes that we are using. Let rj(t) instead be a process with cumulant 
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generating functional W[A(t)] so that the cumulants of r)(t) (which may depend upon x(t)) 
are given by functional derivatives with respect to X(t). This process will have its own action 
S[r](t)] and the path integral can be written as 



, e s[v(t)] 



P[x(t)\y,t ] = Jvr ] (t)S[x{t)-f(x,t)-r ] (t)-y5(t-t )} l 

T>r)(t)Vx(t)e~ S Ht)(.i(t)-f(.^)-yS{t-t ))+mv{t)dt-s[ v (t)] 



Noting that 

' ' Vr](t)e j "*(*)»?(*)#-%(*)] = e w[x(t)] 



is the definition of the cumulant generating functional for rj(t), we have that the path integral 
can be written as 



P[x(t)\y,t ] = J p ?7 (t)Px(t) e -/ 5 ( t )(*(*)-/(*,*)-y*(t-to))*+w[x(t)] 
In the cases where the input rj(t) is delta-correlated in time, we obtain 

oo „ oo „ 

mm = E / 9n(x(t))x(t) n dt = J2 / £n ^ xm ^ dt 

n=l n=l,m=0 

where we have Taylor expanded the functions g n {x). For example, the Ito process above 
gives 

W[x(t)} = °J x{tfdt 
i.e. t>2o = D and all other v nm = 0. 

A. Ornstein-Uhlenbeck Process 

Consider the Ornstein-Uhlenbeck process 

x(t) + ax(t) - VDq{t) = 
with initial condition x(0) = y. The action is 



S[x, x] 



x(t) (x(t) + ax(t) - yS(t - t )) - y£ 2 (t) 



dt 
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Defining an inverse propagator 



G-\t-t')=[j t +a)5{t-t') 



the action is 



S[x, x] = J x(t)G~ l (t - t')x(t')dtdt' - J yx(t)S(t - t )dt - J ®x(t) 2 dt 

and the generating functional is 

Z[J,J] = J p ar ( t )2?x(t) e -^^+/J(0*W*+/J(*)*(*)«e 

This path integral can be evaluated directly as a Gaussian integral since the action is 
quadratic. In fact integrating by x(t) gives the Onsager-Machlup path integral [21 EI], which 
will have a Jacobian factor depending upon whether we use Ito, Stratonovich, or some other 
convention for our SDE. With the Ito convention, this Jacobian is 1. However, the gener- 
ating functional can also be evaluated by expanding the exponent around the "free" action 
given by Sp[x(t),x(t)] = f x(t)G _1 (t — t')x{t')dtdt' . We will demonstrate this method since 
it forms the basis for perturbation theory for non-quadratic actions. Expand the integrand 
of the generating functional as 

Z[J, J} = J Vx(t)Vx(t)e-!^*®^W*M (l + ii + + ^/i 3 + • • •) (7) 

where 

H = y j x(t)5{t-t )dt + J ®x 2 (t)dt + J J(t)x(t)dt + J J(t)x{t)dt 

The generating functional is now expressed as a sum of moments of the free action, which 
are calculated from the free generating functional 

Z F [J,J] = J Vx{t)Vx{t)e~ ^ dtdt ' i{t)G ' 1{t ' t ' )x{t ' )+ 

From equation ([2]), we can write this as 

Z F [J i 3\ = ef i ® G W™ (8) 
where G(t,t') is the operator inverse of i.e. 

j dt"G-\t,t")G(t",t') = + aj G(t,t') = S(t - t') 
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Therefore 

G(t,t') = H(t - lf)e- ait ~ e) 

where H(t) is the left continuous Heaviside step function (i.e. H(0) = 0, lim t _ s>0 + H(t) = 1 
and thus lim ti _ >t + G(ti,t 2 ) = 1, G(t,t) = 0). The choice of i?(0) = is consistent with the 
Ito condition for the SDE and insures that the configuration variable x(t) is uncorrelated 
with future values of the stochastic driving term. Other choices for H(0) represent other 
forms of stochastic calculus (e.g. H(0) = 1/2 is the choice consistent with Stratonovich 
calculus) [30]. The free moments are given by 

since Zp[0, 0] = 1. We use a subscript F to denote expectation values with respect to the 
free action. From the action of ([8]), it is clear the nonzero free moments must have equal 
numbers of x(t) and x(t) due to Wick's theorem, which applies here for contractions between 
x(t) and x(t). For example, one of the fourth moments is given by 

(z(*i)z(* 2 )x(*3)x(*4)>f = G(t h t 3 )G(t 2 ,h) + G(t ll U)G(t 2 ,t 3 ) 

Now the generating functional for the OU process ^ can be evaluated. The 
only surviving terms in the expansion will have equal numbers of x(t) and x(t). 
Thus only terms with factors of J x(to)J(ti)x(ti)dti, (D/2) f x 2 {ti)J 2 {t2)x 2 {t2)dtidt2 and 
J J(ti)x(ti)J(t2)x(t2)dtidt2 (and combinations of the three) will survive. For the OU pro- 
cess, the entire series is summable. First consider the case where D = 0. Because there must 
be equal numbers of x(t) and x(t) factors in any non-zero moment due to Wick's theorem, 
in this case the generating functional has the form 



Z 



H ~n / ( II J(tiMtiHtj)[y$(tj - to) + J(tj)] ) f[ dudt, (9) 

m=\ \i,j=l I p ij=l 



From Wick's theorem, the free expectation value in Q will be a sum over all possible 
contractions between x(t) and x(t) leading to ml combinations. Thus ^ is 



y- 

^ ml 

m=l 



y J J(t 1 )G(t 1 ,t )dt 1 + J J(t')J(t")G(t',t")dt'dt" 
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which means the series is an exponential function. The other term in the exponent of (10) 
can be similarly calculated resulting in 

Z[J(t),J(t)]=exp (y J J(t 1 )G(t 1 ,t )dt 1 + J J^) J(i 2 )G(t 1? ^dt^ 

+f J J(t 1 )J(t 2 )G(t 1 ,t")G(t 2 ,t")dt"dt 1 dt 2 ^j (10) 

The cumulant generating functional is 

W[J(t),J(t)] = y J J(t)G{t,t )dt+ J J(t')J(t")G(t',t")dt'dt" 

+ Y / J(t')At")G(t' ,t)G(t" ,t)dtdt'dt" (11) 

The only nonzero cumulants are the mean 

(x(t)) =yG{t,t ) 

the response function 

5 5 

(x(h)x(t 2 )) c = 6J{h) W[J,J] J= j =0 = G(t 1 ,t 2 ) 



and covariance 



{x{h)x{t 2 )) c = (x(h)x(t 2 )) - (x(ti))(x(t 2 )> 
6 6 

= 6J(t 1 )5J(t 2 ) W[J,J]j=J=0 
= D G{t x ,t)G{t 2 ,t)dt 



Closed form expressions for the cumulants are obtained by using the solution for the 
propagator G. Hence, the mean is 

(x(t)) = ye- a{t - to) H(t - t ) (12) 

the response function is 

(x{h)x(t 2 )) = e'^-^Hih - t 2 ) 

and the covariance is 

(x(ti)x(t 2 )) c = D f 2 e-^-^e-^-^Hih - t')H(t 2 - t')dt' 

J t 



14 



For t 2 > ti > t 

e 2a(ti-t 2 ) _ e -a(*i+t 2 -2t ) 

(x(tMt 2 )) c = D 

For tx = t 2 = t 



{<tf)c = £(1 - e~**) (13) 



IV. PERTURB ATIVE METHODS AND FEYNMAN DIAGRAMS 

If the SDE is nonlinear, then the generating functional cannot be computed exactly as 
in the linear case. However, propagators and moments can be computed perturbatively. 
The method we use is an infinite dimensional generalization of Laplace's method for finite 
dimensional integrals [21] . In fact, the method was used to compute the generating functional 
for the Ornstein-Uhlenbeck process. The only difference is that for nonlinear SDEs the 
resulting asymptotic series is not generally summable. 

The strategy is again to split the action £[2;, x] = Sp + Sj, where Sp is called the "free" 
action and Sj is called the "interacting" action. The generating functional is 

^JH/^Ke-sMWW- (!4) 

The moments satisfy 



libit \ _, 1 1 it 1 b 

nn^)^))=^nn 



m n r- r 





and the cumulants satisfy 

m n r. j- 

3 ' C ^ 



(15) 



j=j=o 



(16) 



j=j=o 



The generating functional is computed perturbatively by expanding the integrand of (14) 
around the free action 

Z[J, J] = J VxVie- SF[x ' i] ^1 + 5/ + J Jxdt + J Jxdt 

j J~xdt^j +is / 3 + ---j 



+^7 ( SV + / Jxdt + 



2! 

Hence, the generating functional can be expressed in terms of a series of free moments. 
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Now we apply this idea to the example nonlinear SDE 

x = —ax + bx 2 + yS(t — t ) + \fDx^r]{t) 

for some n > 0. For example, n = corresponds to standard additive noise (as in the OU 
process), while n = 1 gives multiplicative noise with variance proportional to x. The action 
for this equation is 

D 



2 



= S F [x,x] - yx(t ) -b J dtx(t)x 2 (t) - J dtx 2 x n ^ (17) 

where we have implicitly defined the "free" action Sp[x,x] = f dtx(x + ax). Expectations 
with respect to this free action are 

(x(t)x(t')) F = G(t,t') 

where the propagator obeys 

(| + «) G(t,t') = 6(t-t') 
and all other moments are zero. The generating functional is 

The Taylor expansion of the exponential around the free action gives 
Z[J, 1} = / D^-'M (> • ,/ ^ * + iM y + f / *V * + / Jxdt + / 



^ (b J xx 2 dt + x(t )y +j J x 2 * 11 dt + J Jxdt + J Jxdt^j 



2 



+ — | b I .{■.(- r/M .(■(/„)(/ I .r.rdl \- j J.r<tt+ I J.rdt ) 

Because the free action Sf is bilinear in x, x, the only surviving terms in the expansion are 
those with equal numbers of x and x factors. Also, because of the Ito condition, H(0) = 0, 
these pairings must come from different terms in the expansion, e.g. the only term surviving 
from the first line is the very first term, regardless of the value of n. All other terms come 
from the quadratic and higher terms in the expansion. For simplicity in the remainder of 
this example we limit ourselves to n — 0. Hence, the expansion includes terms of the form 
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i 3\ 72 d [.., r _ 2 , /• _ 2 , 

+ — — o — I xx at I xx at I x at 

+ 1 1 § / i2 dt I Jxdt J J * dt + h™f I ix ' dt I i2dt I ~ Mi 

+ ^7yb J xx 2 dt(x(to)y) 2 J Jxdt + ^—(x(t )yf j Jxdt J Jxdt 
+ — 5!£»-^- J xx 2 dt J x 2 dtx(t )y J Jxdt J Jxdt + - 

Note that this not an exhaustive list of terms up to fifth order. Many of these terms will 
vanish because of G(t,t') oc H(t — t') and H(0) = 0. The combinatorial factors arise from 
the multiple ways of combining terms in the expansion. There are n! ways of combining 
terms at order n and terms with m repeats are divided by a factor of m\. Completing the 
Gaussian integrals using Wick's theorem then yields 

Z[J,J] = Z F [0,0](1 
+ y G(t u t )J(t 1 )dt 1 + / J(t 1 )G(t h h)J(t 2 )dt 1 dt 2 



+ D j GfatjGfatJJ^J^dhdhdts 
+ bD J G{t^t 2 ) 2 G{h,t{)J{t z )dt^dt2dt z 
+ by 2 J G(ti,to) 2 G(t 2) ti)J(t2)dtiG?t2 
+ y 2 J GitutoWtJdt! J G{t 2 ,t )J{t 2 )dt 2 

+ 2bDy J G{h, t 2 )G(ti, t )G{t 3 , h)G{U, t 2 )J(t 3 )J{t 4 ) dt x dt 2 dt z dt± + ■ ■ ■) (18) 
As above, we have Z F [0, 0] = 1. 



The moments and cumulants are obtained from (15) and (16) respectively. For example, 
the mean is given by 

(x(t)} = * . - Z[J(t), J(t)\ 7 _ n f _ a 
V V " Z[0,0] SJ(t) )\J-o,J-o 

= yG(t, t ) + bD J G(t, ti)G(ti, t 2 ) 2 dt x dt 2 + by 2 J G(t, ti)G(*i, t ) 2 dh + ■ ■ ■ (19) 
The covariance is 

S S 

(x(s)x(t)) = ^Z[J(t), J(t)] 7 _ n ?_ n 

XV ; K " 5J(s)5J{t) 1 W V ;jJ -°' J -° 
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= D J G(s,t 1 )G(t,t 1 )dt 1 + y 2 G(s,t )G(t,t ) 

+ 2bDy J G(t 1 ,t 2 )G(t 1 ,t )G(s,t 1 )G(t,t 2 )dt 1 dt 2 

+ 26L>y y G(t 1 ,t 2 )G{t 1 ,to)G(t,t 1 )G(s,t 2 )dt 1 dt2--- 

The first cumulant is the same as the mean but the second cumulant or covariance is 

S S 

(x(s)x(t)) c = - - \nZ[J(t),J(t)] J=oJ=0 
dJ(s) dJ(t) 



Z 6J(s) S J{t) 



5 5 
-Z^—Z 



=o,j=o 5J W 



J=0,J=0 



= D y G(s,*i)G(Mi)<**i 

+ 2bDy y G(t 1 ,t 2 )G(t 1 ,t )G(s,t 1 )G(t,t 2 )dt 1 dt2 

+ 26^ y G(ti,t 2 )G(t 1 ,t )G(t, t^Gis, t 2 ) dt x dt 2 ■ ■ ■ (20) 

As can be seen in this example, the terms in the perturbation series become rapidly 
unwieldy. However, a convenient means to keep track of the terms is to use Feynman 
diagrams, which are graphs with edges connected by vertices that represents each term in 
the expansion of a moment. The edges and vertices represent terms (i.e. interactions) in 
the action and hence SDE, which are combined according to a set of rules that reproduces 
the perturbation expansion shown above. These are directed graphs (unlike the Feynman 
diagrams usually used for equilibrium statistical mechanics or particle physics). The flow 
of each graph, which represents the flow of time, is directed from right to left, points to 
the left being considered to be at times after points to the right. The vertices represent 
points in time and separate into two groups: endpoint vertices and interior vertices. The 
moment (^i = ix{U) YljLi^i^j)^ * s represented by diagrams with iV final endpoint vertices 
which represent the times ti and M initial endpoint vertices which represent the times tj. 
Interior vertices are determined from terms in the action. 

Consider the interacting action expressed as the power series 

,.nm poo 

Si= v ™= E —r dt ~ xnxm ( 21 ) 

n>2,m>0 n>2,m>0 ' 

where n and m cannot both be < 1 (those terms are part of the free action). (Non- 
polynomial functions in the action are expanded in a Taylor series to obtain this form.) 
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There is a vertex type associated with each V nm . The moment (IIi=i x (ti) Ylj=i x (tj)) 
is given by a perturbative expansion of free action moments that are proportional to 
X {U) YljLi x{tj)V(N v )^ where V(N V ) represents a product of N v vertices. Each 
term in this expansion corresponds to a graph with N v interior vertices. We label the fcth 



vertex with time As indicated in equation (21), there is an integration over each such 
interior time point, over the interval (io,oo). The interaction V nm produces vertices with 
n edges to the left of the vertex (towards increasing time) and m edges to the right of 
the vertex (towards decreasing times). Edges between vertices are due to Wick's theorem, 
which tells us that every x(t') must be joined by a factor of x(t) in the future, i.e. t > t', 
because G(t,t') oc H(t — t'). Also, by the Ito condition H(0) = 0, each edge must connect 
two different vertices. All edges must be connected, a vertex for the interaction V nm must 
connect to n edges on the left and m edges on the right. 

Hence, terms at the N v th order of the expansion for the moment (n^i X {U) Y\j = \ %(tj 



are given by directed Feynman graphs with N final endpoint vertices, M initial endpoint 
vertices, and N v interior vertices with edges joining all vertices in all possible ways. The 
sum of the terms associated with these graphs is the value of the moment to N v th order. 
Figure [l] shows the vertices applicable to action (17) with n — 0. Arrows indicate the flow 
of time, from right to left. These components are combined into diagrams for the respective 
moments. Figure [2] shows three diagrams in the sum for the mean and second moment 
of x(t). The entire expansion for any given moment can be expressed by constructing the 
Feynman diagrams for each term. 

The application of Feynman diagrams for computing the diagrams corresponding to terms 
in a perturbative expansion is encapsulated in the Feynman rules: 

A) For each vertex type, V nm , which appears k times in the diagram, there is a factor of 



fc! • 



B) For n distinct ways of connecting edges to vertices that yields the same diagram, i.e. 
the same topology, there is an overall factor of n. This is the combinatoric factor from the 
number of different Wick contractions that yield the same diagram. 

C) Each vertex interaction V nm adds a factor — ^r- The minus sign enters because the 
action appears in the path integral with a minus sign. 

D) For each edge between times t and t', there is a factor of G(t,t'). 
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a) b) 




c) d) 




FIG. 1: Feynman diagram components for a) an edge, the propagator G(t,t'), and vertices b) 
J bxx 2 dt, c) f yx5(t — to)dt, and d) J ^x 2 dt. 

E) There is an integration over the times t of each interior vertex over the domain (to, oo). 
Comparing these rules with the diagrams in Figure [2j one can see the terms in the 



expansions in equations ( 19 ) and (20 ), with the exception of the middle diagram in Figure|2p. 
An examination of Figure [2^, shows that this middle diagram is two copies of the first 
diagram of the mean. Topologically, the diagrams have two forms. There are connected 
graphs and disconnected graphs. The disconnected graphs represent terms that can be 
completely factored into a product of moments of lower order (cf. the middle diagram in 
Figure^). Cumulants consist only of connected graphs since the products of lower ordered 
moments are subtracted from the moment by definition. Thus, moments and cumulants can 
be computed directly from the diagrams that represent them. The connected diagrams in 
Figure [2] lead to the expressions (19) and (20). In the expansion (18), the terms that do not 



include the source factors J and J only contribute to the normalization Z[0, 0] and do not 



affect moments because of (15). Borrowing terminology from quantum theory, these terms 
are called vacuum graphs and consist of closed graphs, i.e. they have no initial or trailing 
edges. In the cases we consider, all of these terms are 0, which implies Z[0, 0] = 1. 

The diagrammatic expansion is particularly useful if the series can be truncated so that 
only a few diagrams need to be computed. There are two types of expansions depending 
on whether the nonlinearity is small or the noise source is small. In quantum theory, the 
small nonlinearity expansion is called a weak coupling expansion and the small fluctua- 
tion expansion is called a semiclassical or loop expansion. The weak coupling expansion is 
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a) 



b) 




FIG. 2: Feynman diagrams for a) the mean and b) second moment. 

straightforward. Suppose one or more of the vertices is associated with a small parameter a. 
These vertices define the interacting action Si as demonstrated above. Each appearance of 
that particular vertex diagram contributes a factor of a and the expansion can be continued 
to any order in a. 

For the loop expansion, let us introduce the factor h into the generating functional: 

Z[J, J] = J Pa;(t)I>5(t)e-K flr [ a! W^W]-/-'W*(*)*-/J(t)*(*)«a) (22) 

According to the Feynman rules described above, with this change each diagram gains a 
factor of h for each edge (internal or external) and a factor ofl/h for each vertex. Let E be 
the number of external edges, I the number of internal edges, and V the number of vertices. 
Then each connected graph now has a factor h I+E ~ v . It can be shown via induction that 
the number of closed loops L in a given connected graph must satisfy L = I — V + 1 [19J. To 
see this note that for diagrams without loops any two vertices must be connected by at most 
one internal edge. Since the diagrams are connected we must have V — I + 1 when L = 0. 
Adding an internal edge between any two vertices increases the number of loops by precisely 
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one. Thus we see that the total factor for each diagram may be written h E+L ~ 1 . We can 
organize the diagrammatic expansion in terms of the number of loops in the graphs. This 
is called the loop expansion. For the mean which has one external edge there are no factors 
of h at lowest order. Higher cumulants (which are determined by connected graphs) gain 
additional factors of h, e.g. the variance goes as h at lowest order in the loop expansion. 

Loop diagrams arise because of nonlinearities in the SDE that couple to moments of the 
driving noise source. For example, the middle graph in Figure [2^, describes the coupling 
of the variance to the mean through the nonlinear x 2 term. This produces a single loop 
diagram which is of order h, compared to the order 1 "tree" level mean graph. Compare 
this factor of h to that from the tree level diagram for the variance, which is order h. This 
same construction holds for higher nonlinearities and higher moments for general theories. 
The loop expansion is thus a series organized around the magnitude of the coupling of higher 
moments to lower moments. 

As an example, consider the action 

S[x,x} — J x(x — f(x(t),t)) — a 2 -x 2 g 2 (x(t),t) dt 

where a is a small parameter and / and g are of order one. Rescale the action with the 
transformation x — )■ xja 2 and J — > J /a 2 . The rescaled action now has the form 

S[x, x] = ^ Jx(x- f(x(t),t)) - h 2 g 2 (x(t),t) dt 

The generating functional is 

Z[J : J\ = [ p x p§ e -^(/^( :i -/( £t; (*).*))-| £ V(a:(t),t)'it+/j : xdt+/Jxdt) 



The loop expansion in this construction is explicitly a small noise expansion because a 2 
plays the role of h in the loop expansion. 

Now consider the one loop correction to the linear response, (x(t)x(t')) , when n — 1 in 



action (17). For simplicity, we will assume the initial condition y = 0. In this case, the 
vertex in Figure [T|i now appears as in Figure [3] The linear response (x(t)x(t')) will be given 
by the sum of all diagrams with one entering edge and one exiting edge. At tree level, there 
is only one such graph, equal to G(t,t'). At one loop order, we can combine the vertices in 
Figures [TJo and [l]i to get the second graph shown in Figure [4] to obtain 

(x{t)x{t')) = G(t,t') + bD J dt 1 dt 2 G( y t,t2)G(t2,t 1 ) 2 G( y t 2 ,t') 
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FIG. 3: Vertex for multiplicative noise with n = 1 in the action (17). This vertex replaces the one 
in Figure fTH. 



+ 




FIG. 4: Feynam diagrams for the linear response, (x(t)x(tf)), to one loop order. 



t-t' 



-a(t-t') 



This loop correction arises because of two types of vertices. There are vertices that we 
call "branching" (as in Figure [3]), which have more exiting edges then entering edges. The 
opposite case occurs for those vertices which we call "aggregating" . Noise terms in the SDE 
produce vertices with more than one exiting edge. As can be seen from the structure of the 
Feynman diagrams, all moments can be computed exactly when the deterministic part of 
the SDE is linear because it only involves convolving the propagator (i.e. Green's function) 
of the deterministic part of the SDE with the driving noise term, as in the case of the OU 
process above. On the other hand, nonlinearities give rise to vertices with more than one 
entering edge. 

If there are no branching vertices in the action (i.e. terms quadratic or higher in x), we 
do not even have an SDE at all, but just an ordinary differential equation. Consider the 



expansion of the mean for action (17) for the case where D = (so that there is no noise 
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term). From equation (19), we have 

(x{t)) = yG(t,t ) + by 2 I G{t, t f dt x + 



The expansion for D = will be the sum of all tree level diagrams. It is easy to see that 
in general this expansion will be the perturbative expansion for the solution of the ordinary 
differential equation obtained by discarding the stochastic driving term. In other words, the 
sum of the tree level diagrams for the mean satisfies 

j t (x(t)) tiee = -a(x(t)) tree + b(x(t)) 2 tiee (23) 

along with the initial condition (x(£o))tree = V- Similarly, the sum of the tree level diagrams 



for the linear response, (x(t)x(t')) tTee = G tree (t } t') } is the solution of the linearization of (23) 
with a Dirac delta functional initial condition, i.e. the propagator 

^-G tIee {t, = -aG tiee (t, t') + 2b(x{t)) tiee G tTee {t, t') + 5(t - if) 
at 

The semiclassical approximation amounts to a small noise perturbation around the solution 
to this equation. We can represent the sum of the tree level diagrams graphically by using 
bold edges, which we call "classical" edges, as in Figure |5j We can then use the classical 
edges within the loop expansion to compute semiclassical approximations to the moments 
of the solution to the SDE. The one loop semiclassical approximation of the mean for the 
case n = is given by the sum of the first two graphs in Figure |2^i with the thin edges 
replaced by bold edges. For the covariance, the first graph in Figure [2]d suffices, again with 
thin edges replaced by bold edges. These graphs are equivalent to the equations: 



(x(t)} = (x(t)) tIee + bD / dh / dt 2 G tIee (t,t 2 )G tIee (t2,tiY (24) 

J to J to 

and 

rmm(t,t') 

(X{t)x{t')) = D dtlGtrecfotOGtreet*',*!) ( 2 5) 

Jtn 



Using equation ( 25 ) in ( 24 ) gives 



(x{t)) = (x{t)) tree + bD [ dt 1 G trec {t,t 2 )(x{t 2 )x{t 2 )) 
This approximation is first order in the dummy loop parameter h for the mean (one loop) 



and covariance (tree level). For the case n — 1, equations (24) and (25) are 

(x(t)) = (x(t))trco + bD i dt x / dt 2 G tTce {t,t 2 )G trcc {t 2 ,ti) 2 (x(t)) tr . ec 

J t J t 
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= + ^ + ^ + ^ ^ — + ... 

FIG. 5: Bold edges represent the sum of all tree level diagrams contributing to that moment. Top) 
the mean (x(t)) t rec- Bottom) linear response Gtree(t,t'). 

and 

rmm(t,t') 

(x{t)x{t')) = D G , tree (t,f 1 )G tree (t / ,t 1 )(a;(i)) tree 

J t 

Using the definition of Gr t ree(£> the self-consistent semiclassical approximation for (x(t)) 
to one-loop order is 

d r* 

-(x(t)} + a(x(t)} - b(x(t)} 2 = bD j dtiG trcc (Mi) 2 (x(t)) 

or 

d r* 

-{x(t)} + a(x(t)} - b{x(t)} 2 = b j dhixfaWh^c 

The semiclassical approximation known as the "linear noise" approximation takes the tree 
level computation for the mean and covariance. The formal way of deriving these self- 
consistent equations is via the effective action, which is beyond the scope of this review. We 
refer the interested reader to [TU] . 



V. CONNECTION TO FOKKER-PLANCK EQUATION 

In stochastic systems, one is often interested in the PDF p(x, t), which gives the probabil- 
ity density of position x at time t. This is in contrast with the probability density functional 
P[x(t)] which is the probability density of all possible functions or paths x(t). Previous 
sections have been devoted to computing the moments of P[x(t)}, which provide the mo- 
ments of p(x,t) as well. In this section we leverage knowledge of the moments of p(x, t) 
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to determine an equation it must satisfy. In simple cases, this equation is a Fokker-Planck 
equation for p(x, t). 

The PDF p(x, t) can be formally obtained from P[x(t)} by marginalizing over the interior 
points of the function x(t). Consider the transition probability U(xi,ti\xo,to) between two 
points xo, to and xi, ti. This is equal to p(x, t) given the initial condition p(x, t ) = 5(x — x ). 
In terms of path integrals this can be expressed as 

"(x(ti)=xi) 



P(X(tl)=Xl) 

U^t^to) = / Vx{t)P[x{t)\ 

J (x(ta)=xn) 



'(x(t )=x ) 

where the upper limit in the integral is fixed at x{t\) = x\ and the lower at x(to) = x . The 
lower limit appears as the initial condition term in the action and can thus be considered 
part of P[x(t)}. The upper limit on the path integral can be imposed with a functional Dirac 
delta via 

U(x 1 ,t 1 \x ,t ) = J Vx(t)5(x(t 1 )-x 1 )P[x(t)] 
which in the Fourier representation is given by 

U(x 1 ,t 1 \x , t ) = ^-jd\j Vx(t) e x(x ^-^P[x(t)] 

where the contour for the A integral runs along the imaginary axis. This can be rewritten 

as 

^o, *o) = ^~j d\e- x ^- Xo) Z CM (\) (26) 
in terms of an initial condition centered moment generating function 



Z CM (\) = J Vxe x ^- Xo) P[x{t)\ 



where the measure Vx(t) is defined such that Z C m(0) = 1. Note that this generating func- 
tion Z CM (\) is different from the generating functionals we presented in previous sections. 
-Zcm(A) generates moments of the deviations of x(t) from the initial value x at a specific 
point in time t. Taylor expanding the exponential gives 



Zcu(X) = l + J2 -i An <(^0 - x oT)x(t )= Xo 
1 

n=l 

where 

- x ) n ) x{ t o)=xo = J Vx (x(*i) - x ) n P[x(t)\ 
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Inserting into (26) gives 



Using the identity 



i/ d Ae«A»=(-A)"^ 1 - Io) 



{(x(ti) - x ) n ) x (t )= X0 S(xi - x ) (27) 



results in 

too 1 / d \ n 
n=l x ' 

The probability density function p(y, t) obeys 

p{y, t + At) = J U{x, t + At\y', t)p{y', t)dy' (28) 



Inserting Q2TJ ) gives 

/ °° 1 / r) \ n 

" ' ((x^+Atj-i/n^j^f) 



dy, 



Expanding p(y, t + At) and the moments in a Taylor series in At gives 

«p(v, t} At = (-^ " 1 ((*(* + At) - y rUt)= y p(y, t) + 0(At 2 ) 



at 

71=1 



since x(t) = y. In the limit At — > we obtain the Kramers-Moyal expansion 
dp{y,t) f 1 / 9 



at 

ra=l 



1 / f) \ n 

^\-dy) D n(y,t) P (y,t) + o(At 2 ) 



where the jump moments are defined by 

((x(t + At) - y) n ) 



D n (y,t) = lim 



(29) 

At -° A * x(t)=v 

As long as these limits are convergent, then it is relatively easy to see that only connected 
Feynman graphs will contribute to the jump moments. In addition, we can define z = x —y, 
z = x and use the action S[z(t) + y, z(t)]. This shift in x removes the initial condition term. 
This means we can calculate the nth jump moment by using this shifted action to compute 
the sum of all graphs with no initial edges and n final edges (as in Figure [TJd for n = 2). 
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As an example, consider the Ito SDE ([3]). From the discretization (J5J), where h = At, it 
is found that 



At->0 At 



fi(y)At - gi(y)wiVAt 
lim /_X (30) 



Which yields £>i(y,t) = f(y,t), D 2 = g{y,tf and D n = for n > 2. Thus for the Ito SDE 
(j3|, the Kramers- Moyal expansion becomes the the Fokker-Planck equation 
dp(y,t) ( d „ , . lfl 3 n . A , . 



<9i \ dy ' 2 (9 2 y 
We have D n = for n > 2 even though there are non-zero contributions from connected 
graphs to these moments for n > 2 in general. However, all of these moments require the 
repeated use of the vertex with two exiting edges; this will cause D n oc At m for some m > 1 
and thus the jump moment will be zero in the limit. 

We can envision actions for more general stochastic processes by considering vertices 
which have more than two exiting edges, i.e. we can add a term to the action of the form 

SV^m^] = — J dtx n h{x) 

for some n and function h(x). This will produce a non-zero D n . The PDF for this kind of 
process will not in general be describable by a Fokker-Planck equation, but will need the full 
Kramers-Moyal expansion. If we wished to provide an initial distribution for x(to) instead 
of specifying a single point, we could likewise add similar terms to the action. In fact, the 
completely general initial condition term is given by 

Siaitiai[£(<o)] = - lnZ y [x(t )] 

where Z y is the generating functional for the initial distribution. In other words, the initial 
state terms in the action are the cumulants of the initial distribution multiplied by the 
corresponding powers of x(t ). 

Returning to the Ito process ([3]), the solution to the Fokker-Planck equation can be 



obtained directly from the path integral formula for the transition probability (26). Let 
lnZ[A] be the cumulant generating function for the moments of x(t) at time t. It can be 
expanded as 



^cm[A] = exp 



n=l 
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yielding 



p(x,t) 



1 



dXe Xx exp 



£^ B <*(*) B > 

n=l 



C 



For the Ornstein-Uhlenbeck process the first two cumulents are given in (12) and (13) 
yielding (assuming initial condition x(to) = y) 



a f _ a ( x _ ye -q(t-* ))2 ^ 

?WJ - \l 7rD ( 1 _ e -2a(t-t )) eXP ^ £)(1 _ e -2o(t-to) I 



VI. FURTHER READING 



The reader interested in this approach is encouraged to explore the extensive literature on 
path integrals and field theory. The reader should be aware that most of the references listed 
will concentrate on applications and formulations appropriate for equilibrium statistical 
mechanics and particle physics, which means that they will not explicitly discuss the response 
function approach we have demonstrated here. For application driven examinations of path 
integration there is Schulman[25] and Kleinert [H] . More mathematically rigorous treatments 
can be found in Simon [26] and Glimm and Jaffe [27]. For the reader seeking more familiarity 
with concepts of stochastic calculus such as Ito or Stratonovich integration there are applied 
approaches [3] and rigorous treatments [28] as well. Zinn- Justin [19] covers a wide array 
of topics of interest in quantum field theory from statistical mechanics to particle physics. 
Despite the exceptionally terse and dense presentation, the elementary material in this 
volume is recommended to those new to the concept of path integrals. Note that Zinn- Justin 
covers SDEs in a somewhat different manner than that presented here (the Onsager-Machlup 
integral is derived; although see chapters 16 and 17), as does Kleinert. We should also point 
out the parallel between the form of the action for exponential decay (i.e. D = in the 
OU process) and the holomorphic representation of the harmonic oscillator presented in 
[19] . The response function formalism was introduced by Martin, Siggia, Rose [22]. Closely 
related path integral formalisms have been introduced via the work of Doi |6] and Peliti 
[7J which have been used in the analysis of reaction-diffusion system [8 HTCTI 123] . Uses of path 
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integrals in neuroscience have appeared in [T3HT6] . 
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